! VSemiletov, 08.2011
subroutine AcousticSurfInit()
use Acoustics, only: LeftDiscActive,RightDiscActive, CenterDir, Length, factor, &
                     nAcousticSurfaces,AcousticSurface, &
                     Sphere, Center, Radius
use Data_cfoam, only: meshCellOut,nCellsOut,areaOut,xOut,normalOut,nUnusedSections, &
                      processor, ie, nump,rank,  ierr, &
                      nUnusedSections,nUnusedSectionsSum
use Data_Mesh, only: nCells, meshPoints, meshCellPoints, meshFaceCenters, boundaryMeshStart, boundaryMeshSize, &
                     meshFaceOwner, meshSf
implicit none
    include 'mpif.h'

    integer :: i, jj, iCell, j, k, iFace
    integer :: LengthName
    real(8) :: CoorCell(3,8)
    character(50) :: name
    integer :: nCellsOutSum, bmeshs

    nCellsOutSum = 0

    if(nAcousticSurfaces.le.0)then
        if(rank.eq.0) then
            write(*,*) '@Warning: no Acoustics surfaces!'
        end if
        return
    end if

    do i = 1,nAcousticSurfaces
        if(AcousticSurface(i)%DefCellsFromContour)then

            if(AcousticSurface(i)%Sphere)then
                Sphere = .true.
                Center(1:3) = AcousticSurface(i)%Center(1:3)
                Radius = AcousticSurface(i)%Radius
            else
                LeftDiscActive = AcousticSurface(i)%LeftDiscActive
                RightDiscActive = AcousticSurface(i)%RightDiscActive
                CenterDir(1:3) = AcousticSurface(i)%CenterDir(1:3)
                Length = AcousticSurface(i)%Length
                factor = AcousticSurface(i)%factor
                LengthName = 0
                do j = 1,280
                    if(AcousticSurface(i)%fileContourIn(j:j).eq.' ')then
                        LengthName = j
                        exit
                    end if
                end do

                call ReadContourData(AcousticSurface(i)%fileContourIn(1:LengthName))
            end if
            
            meshCellOut = 0
            nCellsOut = 0
    	    areaOut = 0.0d0
    	    xOut = 0.0d0
    	    normalOut = 0.0d0
    	                              
            do iCell=1,nCells
                do jj=1,8
                    CoorCell(1:3,jj)=meshPoints(1:3,meshCellPoints(jj,iCell))
                end do
                call meshCellOut_defn(CoorCell,meshCellOut,nCellsOut,iCell, &
    	                              areaOut(nCellsOut+1), xOut(1:3,nCellsOut+1), &
    	                              normalOut(1:3,nCellsOut+1), nUnusedSections)
            end do

            AcousticSurface(i)%nCellsOut = nCellsOut
            allocate(AcousticSurface(i)%CellsOut(nCellsOut))
            allocate(AcousticSurface(i)%Area(nCellsOut))
            allocate(AcousticSurface(i)%xOut(3,nCellsOut))
            allocate(AcousticSurface(i)%normal(3,nCellsOut))
            
            do j = 1,nCellsOut
                AcousticSurface(i)%CellsOut(j) = meshCellOut(j)
                AcousticSurface(i)%Area(j) = areaOut(j)
                AcousticSurface(i)%xOut(1:3,j) = xOut(1:3,j)
                AcousticSurface(i)%normal(1:3,j) = normalOut(1:3,j)
            end do

            ! AcousticsOut = On
            if(AcousticSurface(i)%CalcAcoustics.eqv..false.)then ! AcousticSurface(i)%AcousticsOut(1:1).ne.' ')then
                call AcousticsMeshOut(AcousticSurface(i)%AcousticsOut,AcousticSurface(i)%LengthAcousticsOut)
            end if

!        else
!            call readCellsOutAcousticSurfs(i,processor(1:ie)//name)
        end if

        if(AcousticSurface(i)%nBoundary.gt.0)then

            j = AcousticSurface(i)%nBoundary
            bmeshs = boundaryMeshStart(j)
!            bmeshf = boundaryMeshStart(j) + boundaryMeshSize(j)
            nCellsOut = boundaryMeshSize(j)
            do k = 1,nCellsOut ! bmeshs,bmeshf
                iFace = k + bmeshs
                meshCellOut(k) = meshFaceOwner(iFace)
                areaOut(k) = sqrt(sum(meshSf(1:3,iFace)**2))
                normalOut(1:3,k) = -meshSf(1:3,iFace) / areaOut(k)
                xOut(1:3,k) = meshFaceCenters(1:3,iFace)
            end do

            AcousticSurface(i)%nCellsOut = nCellsOut
            allocate(AcousticSurface(i)%CellsOut(nCellsOut))
            allocate(AcousticSurface(i)%Area(nCellsOut))
            allocate(AcousticSurface(i)%xOut(3,nCellsOut))
            allocate(AcousticSurface(i)%normal(3,nCellsOut))
            
            do j = 1,nCellsOut
                AcousticSurface(i)%CellsOut(j) = meshCellOut(j)
                AcousticSurface(i)%Area(j) = areaOut(j)
                AcousticSurface(i)%xOut(1:3,j) = xOut(1:3,j)
                AcousticSurface(i)%normal(1:3,j) = normalOut(1:3,j)
            end do

            ! AcousticsOut = On
            if(AcousticSurface(i)%CalcAcoustics.eqv..false.)then ! AcousticSurface(i)%AcousticsOut(1:1).ne.' ')then
                call AcousticsMeshOut(AcousticSurface(i)%AcousticsOut,AcousticSurface(i)%LengthAcousticsOut)
            end if

!        else
!            call readCellsOutAcousticSurfs(i,processor(1:ie)//name)
        end if

        name(1:9) = 'CellSurf_'
        name(10:50) = AcousticSurface(i)%name(1:41)
        call writeCellsOutAcousticSurfs(i,processor(1:ie)//name)

        if (nump.gt.1) then
            AcousticSurface(i)%nCellsOutGlobal = 0
            call MPI_ALLREDUCE(AcousticSurface(i)%nCellsOut, &
                               AcousticSurface(i)%nCellsOutGlobal, &
                               1,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,ierr)
            call MPI_Barrier(MPI_COMM_WORLD,ierr)
        else
            AcousticSurface(i)%nCellsOutGlobal = AcousticSurface(i)%nCellsOut
        end if
        
        if(rank.eq.0)then
            if(AcousticSurface(i)%nCellsOutGlobal.eq.0)then
                write(*,*) '@@ Fatal Error: wrong acoustic surface: nCellsOutGlobal==0! ', AcousticSurface(i)%name
                stop '@@Fatal Error: wrong acoustic surface: nCellsOutGlobal==0!'
            end if
        end if
        
        nCellsOutSum = nCellsOutSum + AcousticSurface(i)%nCellsOutGlobal

        if(AcousticSurface(i)%MeshDataOut)then
            call AcousticsMeshDataOut(AcousticSurface(i)%name, &
                                      AcousticSurface(i)%nameLength, &
                                      AcousticSurface(i)%nCellsOut, &
                                      AcousticSurface(i)%CellsOut, &
                                      AcousticSurface(i)%Area, &
                                      AcousticSurface(i)%xOut, &
                                      AcousticSurface(i)%normal)
        end if
        
        call DefineNextFileName(AcousticSurface(i)%AcousticsOut,AcousticSurface(i)%LengthAcousticsOut,.true.)
    end do

    if (nump.gt.1) then
        nUnusedSectionsSum = 0
        call MPI_ALLREDUCE(nUnusedSections,nUnusedSectionsSum,1,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,ierr)
        call MPI_Barrier(MPI_COMM_WORLD,ierr)
    else
        nUnusedSectionsSum = nUnusedSections
    end if

    if(rank.eq.0)then
        write(*,*) 'Surface: nCellsOut = ',nCellsOutSum, ' nUnusedSections = ',nUnusedSectionsSum
        nUnusedSectionsSum = nUnusedSectionsSum * 100 / nCellsOutSum
        if(nUnusedSectionsSum.ne.0)then
            write(*,*) '@Warning: nUnusedSections = ',nUnusedSectionsSum, '%'
        else
            write(*,*) 'nUnusedSections = 0%'
        end if
    end if

end subroutine AcousticSurfInit

subroutine writeCellsOutAcousticSurfs(iSurf,name)
use Acoustics, only: AcousticSurface
use Data_Mesh, only: meshCellCenters
implicit none
    integer,intent(in) :: iSurf ! number of Acoustic surface
    character(*),intent(in) :: name

    integer :: i
    integer :: OpenFile

    open(11,file=name,IOSTAT=OpenFile)
    if (OpenFile.ne.0) then
	    write(*,*) '@@ Fatal Error (writeCellsOutAcousticSurfs): File open error '//name
	    STOP '@@ Fatal Error!'
    end if

    write(11,'(a50)') AcousticSurface(iSurf)%name
    write(11,*) 'nCellsOut:'
    write(11,*) AcousticSurface(iSurf)%nCellsOut
    write(11,*) 'CellsOut  meshCellCenters(1:3)'

    do i=1,AcousticSurface(iSurf)%nCellsOut
        write(11,'(i15,3f8.3)') AcousticSurface(iSurf)%CellsOut(i), &
                                meshCellCenters(1:3,AcousticSurface(iSurf)%CellsOut(i))
    end do

    write(11,*) 'Area'
    do i=1,AcousticSurface(iSurf)%nCellsOut
        write(11,'(i15,f30.15)') AcousticSurface(iSurf)%CellsOut(i), &
                                 AcousticSurface(iSurf)%Area(i)
    end do

    write(11,*) 'xOut(1:3)'
    do i=1,AcousticSurface(iSurf)%nCellsOut
        write(11,'(i15,3f30.15)') AcousticSurface(iSurf)%CellsOut(i), &
                                  AcousticSurface(iSurf)%xOut(1:3,i)
    end do


    write(11,*) 'normal(1:3)'
    do i=1,AcousticSurface(iSurf)%nCellsOut
        write(11,'(i15,3f30.15)') AcousticSurface(iSurf)%CellsOut(i), &
                                  AcousticSurface(iSurf)%normal(1:3,i)
    end do

    close(11)
    return
end subroutine writeCellsOutAcousticSurfs


subroutine readCellsOutAcousticSurfs(iSurf,name)
use Acoustics, only: AcousticSurface
use Data_Mesh, only: meshCellCenters
implicit none
    integer,intent(in) :: iSurf ! number of Acoustic surface
    character(*),intent(in) :: name

    integer :: i, nCellsOut
    integer :: OpenFile

    open(11,file=name,IOSTAT=OpenFile)
    if (OpenFile.ne.0) then
	    write(*,*) '@@ Fatal Error (readCellsOutAcousticSurfs): File open error '//name
	    STOP '@@ Fatal Error!'
    end if

    read(11,*) ! AcousticSurface(iSurf)%name
    read(11,*) ! 'nCellsOut:'
    read(11,*) AcousticSurface(iSurf)%nCellsOut
    read(11,*) ! 'CellsOut  meshCellCenters(1:3)'

    nCellsOut = AcousticSurface(iSurf)%nCellsOut
    allocate(AcousticSurface(iSurf)%CellsOut(nCellsOut))
    allocate(AcousticSurface(iSurf)%Area(nCellsOut))
    allocate(AcousticSurface(iSurf)%xOut(1:3,nCellsOut))
    allocate(AcousticSurface(iSurf)%normal(1:3,nCellsOut))

    do i=1,AcousticSurface(iSurf)%nCellsOut
        read(11,*) AcousticSurface(iSurf)%CellsOut(i)
    end do

    read(11,*) ! 'Area'
    do i=1,AcousticSurface(iSurf)%nCellsOut
        read(11,*) AcousticSurface(iSurf)%CellsOut(i), &
                   AcousticSurface(iSurf)%Area(i)
    end do

    read(11,*) ! 'xOut(1:3)'
    do i=1,AcousticSurface(iSurf)%nCellsOut
        read(11,*) AcousticSurface(iSurf)%CellsOut(i), &
                   AcousticSurface(iSurf)%xOut(1:3,i)
    end do

    read(11,*) ! 'normal(1:3)'
    do i=1,AcousticSurface(iSurf)%nCellsOut
        read(11,*) AcousticSurface(iSurf)%CellsOut(i), &
                   AcousticSurface(iSurf)%normal(1:3,i)
    end do
    
    close(11)
    return
end subroutine readCellsOutAcousticSurfs


! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine NoiseSugnalInit()
use Acoustics, only: nAcousticSurfaces,AcousticSurface, nListeners,Listener, nNoiseSignals,NoiseSignal
use Data_cfoam, only: processor, ie, nump, ierr !, &
!                      nCells, meshCellOut,nCellsOut,areaOut,xOut,normalOut,nUnusedSections, &
!                      meshPoints, meshCellPoints, &
!                      nUnusedSections,nUnusedSectionsSum
implicit none
!    include 'mpif.h'

    integer :: iSurf, iListener, iNoiseSignal
    integer :: n, nFL

    nNoiseSignals = nAcousticSurfaces * nListeners

    allocate(NoiseSignal(nNoiseSignals))

    do iSurf = 1,nAcousticSurfaces
        do iListener = 1,nListeners
            iNoiseSignal = iListener + (iSurf-1) * nListeners
            NoiseSignal(iNoiseSignal)%iSurface = iSurf
            NoiseSignal(iNoiseSignal)%iListener = iListener

            NoiseSignal(iNoiseSignal)%outputFile(1:3) = 'NS_'
            n = AcousticSurface(iSurf)%nameLength
            NoiseSignal(iNoiseSignal)%outputFile(4:3+n) = AcousticSurface(iSurf)%name(1:n)
            
            nFL = 3+n+1
            NoiseSignal(iNoiseSignal)%outputFile(nFL:nFL) = '_'

            n = Listener(iListener)%nameLength
            NoiseSignal(iNoiseSignal)%outputFile(nFL+1:nFL+n) = Listener(iListener)%name(1:n)

            nFL = nFl + n + 1
            NoiseSignal(iNoiseSignal)%outputFile(nFL:nFL) = '_'

            NoiseSignal(iNoiseSignal)%nameFileLength = nFL + 1
        end do
    end do

!        integer :: nValues=0 ! number of time's points for output
!        real(8) :: Pobs(300000) ! pertubation of pressure


end subroutine NoiseSugnalInit

subroutine AcousticsMeshOut(name,LengthName)
use Data_cfoam, only: processor, ie, &
                      meshCellOut,nCellsOut,areaOut,xOut,normalOut
implicit none
    character(50),intent(in) :: name
    integer,intent(in) :: LengthName

    character(55) :: nameMesh
    integer :: i,j
    integer :: OpenFile
    
    nameMesh(1:LengthName) = name(1:LengthName)
    nameMesh(LengthName+1:LengthName+4) = 'Mesh'

    open(11,file=processor(1:ie)//nameMesh(1:LengthName+4),form='unformatted',IOSTAT=OpenFile)
    if (OpenFile.ne.0) then
	    write(*,*) '@@ Fatal Error (AcousticsMeshOut): File open error ', processor(1:ie)//name
	    STOP '@@ Fatal Error!'
    end if

    write(11) nCellsOut
    write(11) (areaOut(i), i=1,nCellsOut)
    write(11) ((xOut(j,i), j=1,3), i=1,nCellsOut)
    write(11) ((normalOut(j,i), j=1,3), i=1,nCellsOut)
    write(11) ( meshCellOut(i), i=1,nCellsOut)

    close(11)

end subroutine AcousticsMeshOut

subroutine AcousticsMeshDataOut(name,length, nCellsOut,meshCellOut, Area,xOut,normal)
use Data_Mesh, only: nPoints, meshPoints, meshCellPoints
use Data_cfoam, only: rank, processor, ie
implicit none
    character(50),intent(in) :: name
    integer,intent(in) :: length
    integer,intent(in) :: nCellsOut
    integer,intent(in) :: meshCellOut(nCellsOut)
    real(8),intent(in) :: Area(nCellsOut), xOut(3,nCellsOut), normal(3,nCellsOut)

    integer :: i,jj
    character(55) :: fileName

    fileName(1:length) = name(1:length)
    fileName(length+1:length+7) = 'meshhsp'

	if(rank.eq.0)then
    	write(*,*) 'File ', name(1:length), ' Writing ...'
    end if    

	open(1,file=processor(1:ie)//fileName(1:length+7),form='unformatted')

	write(1) nPoints,nCellsOut,((meshPoints(jj,i),jj=1,3),i=1,nPoints) &
                    ,((meshCellPoints(jj,meshCellOut(i)),jj=1,8),i=1,nCellsOut)
    close(1)

    if(rank.eq.0)then
    	write(*,*) 'File ', name(1:length), ' was written successfully!'
    end if
    
    fileName(length+1:length+3) = 'hsp'
    open(1,file=processor(1:ie)//fileName(1:length+3),form='unformatted')
    
    write(1) 0d0,0,nCellsOut
    write(1) Area(1:nCellsOut)
    write(1) xOut(1,1:nCellsOut)
    do jj = 1,3
        write(1) (normal(jj,i),i=1,nCellsOut)
    end do

    close(1)

end subroutine AcousticsMeshDataOut